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Using complementary numerical approaches at high resolution, we study the late-time 
behaviour of an inviscid, incompressible two-dimensional flow on the surface of a sphere. 
Starting from a random initial vorticity field comprised of a small set of intermediate 
wavenumber spherical harmonics, we find that — contrary to the predictions of equi¬ 
librium statistical mechanics — the flow does not evolve into a large-scale steady state. 
Instead, significant unsteadiness persists, characterised by a population of persistent 
small-scale vortices interacting with a large-scale oscillating quadrupolar vorticity field. 
Moreover, the vorticity develops a stepped, staircase distribution, consisting of nearly 
homogeneous regions separated by sharp gradients. The persistence of unsteadiness is 
explained by a simple point vortex model characterising the interactions between the 
four main vortices which emerge. 


1. Introduction 

Two-dimensional flow models have often been used as an idealised description of atmo¬ 
spheric and oceanic fluid dynamics. These models afford the maximum simplicity while 
retaining the key advective nonlinearity of fluid motion and the associated scale cascade. 
However, such models fail to capture the intrinsic three-dimensionality of atmospheric 
and oceanic flows at ‘small’ scales, where the effects of rotation and stratification strongly 
compete (Dritsc hel et a/.f l999). Nevertheless, mathematically, much interest remains in 
understanding how nonlinearity organises fluid motion in this simple context. 

In this paper, we examine the behaviour of a bounded two-dimensional incompressible 
inviscid (dissipationless) flow at long times t. The flow evolves freely from some prescribed 
initial condition, which is subsequently transformed by nonlinearity. No forcing is applied. 

This problem has been the subject of much previous research. Almost all laboratory 
experiments and numerical simulations have indicated a final state consisting of a steady 
large-scale circulation pattern filling the whole domain (see Matthaeus et al. | ([1991 b a); 
Montgomery e£ al\ ( 1992[ |l993[ ); Brands et al. (1997); [Bouchet fe Venaille (2012) for a 
sample of the vast literature). Unsteady final states have also been reported in several 
numerical simulations for doubly-periodic planar domains. Segre and Kida found various 


unsteady final states using initial configurations of constant-vorticity patches (Segre & 
Kid a|1998 ), but the simulations appear not to have gone far enough in time for the proper 
final state to evolve (Yin et a/.|2QQ 3). Morita found a late-time persistent oscillatory state 
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by starting from a perturbed zonal flow (Morita 2011). However both examples of un¬ 
steady final states start from special initial conditions that may not develop broad-band 
turbulence in which energy is widely shared among many wavenumbers. It is commonly 
believed that although various metastable states can persist for a long time due to insuf¬ 
ficient mixing, the flow will settle into a steady final state if it undergoes strong turbulent 
mixing during the evolution (for example see Yin et al (|2003|) ). T heoretically, an expla¬ 


nation for this has been proposed by Miller (1990); Miller et al. I (1992) and Robert & 


Sommeria (1991) 

; Robert (1991), hereafter ‘MRS’. (Reviews may be found in Chavanis 

(2002 

); Eyink & Sreenivasan (2006): 

; Majda & Wang ( 

2006); Chavanis ( 

2009); 

Bouchet 

& Venaille 

(2012) 

; Herbert (2013); Qi & Marston (2014 

).) The MRS theory predicts that 


the nonlinear dynamics, characterised by intense turbulent mixing, ultimately results in a 
stationary large-scale, domain-filling flow. Moreover, this flow is obtained as a statistical 
mechanical equilibrium (assuming random fluctuations of fine-scale vorticity) dependent 
only on the inviscid invariants. Those include energy, circulation, and in general an in¬ 
finite set of ‘Casimirs’ that determine the occupation of vorticity uo of each level within 
the fluid (this is the ‘vorticity measure’). 

Verification of the MRS theory has met with mixed success (Q i fc Marston| 2014). The 
theory itself is difficult to work with when many Casimirs are taken into account, and vari¬ 
ous simplifications have been put forth (see e.g. the review by Bouchet & Venaille (2012)). 
Verification is also problematic when nearly all of the numerical methods used suffer from 
some form of dissipation, artificial or real, not present in the original inviscid system for 
which the theory was developed. A few numerical methods have been developed that 
conserve many invariants, including the Casimirs (Abramov & Majda 2003 


Dubinkina 


& Frank 2010). Note that as numerical simulations inevitably have finite resolution, the 


Casimirs or the vorticity field observed numerically are only the coarse-grained correspon¬ 
dence of those quantities of the inviscid system: the resolved vorticity field corresponds 
to the coarse-grained vorticity a), not a;, and the second- and higher-order Casimirs of the 
resolved vorticity are different from the original fine-grained ones dependent on uo. These 
conservative numerical methods would appear to be ideal, and moreover they would 
overcome the long-standing problem of how unresolved scales impact resolved ones, the 
‘closure problem’. However, an irrefutable feature of two-dimensional turbulence is the 
enstrophy cascade, occurring as thin vorticity filaments stretch and thin to ever finer 
scales. By conservation, this cascade is necessary to build up energy at ever larger scales 
through an ‘inverse cascade’ (Fjprtoft 1953 Dritschel et al. 2009 & refs.), as indeed is 
anticipated in the MRS theory. That is, energy growth at large scales requires a cascade 
of enstrophy to small scales. Since the Casimirs involve integrals of powers of the vortic¬ 
ity distribution, all but the first power (the circulation) cannot be even approximately 
conserved at any finite resolution, after sufficiently long time. Forcing conservation of 
quantities which cannot be conserved can only lead to spurious conclusions. 

Verification of the MRS theory is further problematic because the invariant fine-grained 
Casimirs that the theory depends on cannot be measured from any numerical simulation 
of finite resolution. However in practice the resolved Casimirs are often used to determine 
the values of the fine-grained Casimirs (at most approximately). Notably, the coarse¬ 
grained vorticity field uj predicted by the MRS theory often agrees with that emerging at 
large times in numerical simulations, if the late-time values of the resolved Casimirs are 
applied. Still these resolved Casimirs are not conserved and cannot be predicted from an 
arbitrary initial flow state of a numerical simulation, which limits the predicative power 
of the MRS theory. 

Another fundamental issue concerns ergodicity, namely that ensemble averages over 
equi-probable states with the same conserved quantities equate to time averages of a 
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single state over a sufficiently long period of time. Two-dimensional flows are generally 


not ergodic (Whitaker & Turkington 1994 Bricmont et al. 2001 Bouchet & Corvellec 


2010), though they are assumed to be in the MRS theory. It has been argued that MRS 


theory may still apply to states which exhibit strong turbulent mixing over the course of 
their evolution, e.g. to initially small-scale isotropic vorticity fields. There are many such 
states. But it is unclear whether they satisfy ergodicity even approximately. 

The present paper is not intended to be a critique only of MRS theory; rather, we 
wish to present fresh evidence that even the central claim for the existence of a final 
steady equilibrium vorticity distribution uj is likely to be false. Instead of the much 
studied doubly periodic planar domain, we consider flows on the sphere where statistical 
mechanics has not seen much application. While we agree that one can never carry out 
numerical simulations to infinite time, our results strongly suggest that, at least for flows 


on a sphere, unsteadiness persists for all time. Note that both Segre & Kida (1998) and 


Morita (2011) have focused on the type of final unsteadiness that is a result of special 


initial conditions (or insufficiently developed turbulence), but we argue that unsteadiness 
is generic by using a random initial state that generates fairly strong turbulent mixing: 
the unsteadiness is intrinsic to a system of compact coherent vortices. Our simulations 
are carried out much further in time and at much higher resolutions compared to the 
two studies, and the robustness of our result is verified by comparing two very different 
numerical approaches with fundamentally different forms of numerical dissipation. 

The paper is organised as follows. In the next section, we briefly review the equa¬ 
tions and the associated conserved quantities, then describe two very different numerical 
methods employed. Both methods are used in §3 at many different resolutions, on the 
same initial data, to study the long-time convergence of the results. A few conclusions 
are offered in §4. 


2. Governing equations and numerical methods 

2.1. Two-dimensional flow on a sphere 

We wish to study a freely-evolving incompressible two-dimensional flow, with no forcing 
or dissipation, confined to the surface of a unit sphere centered at the origin. In this case, 
the scalar vorticity uj normal to the surface is materially conserved, 


D uj 

Dt 


duj 

dt 


T u • ^ uj — 0 , 


( 2 . 1 ) 


where w = r x is the incompressible velocity field, r is the position of a point on 
the surface (|r| = 1), and ip is the scalar streamfunction. With V • u = 0, the definition 
uj = r • (V x u) yields Poisson’s equation for ip: 


vV = W. 


These equations were originally derived by Zermelo (1902) and Charney (1949). 


( 2 . 2 ) 


2.2. Constants of the motion 


Inviscid two-dimensional flow is an infinite-order Hamiltonian system with the kinetic 
energy 


E = 



(2.3) 


serving as the Hamiltonian. Above, the second expression is obtained after integration 
by parts, and d Tl is the differential surface area element. The integral is taken over the 
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entire surface of the sphere. Besides E, rotational symmetry about any axis implies that 
the vector angular impulse 


// 


ujr dfl 


(2.4) 


is also conserved. Non-zero L implies a mean rotation with vector angular velocity 3L/8n. 
(Below, we only consider flows with L = 0 to avoid starting with a flow organised at the 
largest possible scale.) Finally, material conservation of vorticity and incompressibility 
imply a generally infinite number of additional constraints, namely 


C n = 


//• 


' dn. 


n = 2, 3, 4, 


(2.5) 


are conserved. (C\ = 0 automatically; no net circulation is possible on a closed surface). 
The C n are called the ‘Casimirs’. Their conservation is a consequence of the fact that the 
fractional area occupied by each vorticity value is conserved. For a continuous vorticity 
distribution, this fractional area is generally infinitesimal, and it is then useful to describe 
it as a ‘measure’, or probability density p{uS): the area occupied by vorticity in the range 
[w,uj + dc u\ is 47rp(cj)do;. The conservation of all C n is equivalent to the conservation of 
p(u). 

The quadratic Casimir C 2 is proportional to the ‘enstrophy’ Z (i.e. Z = C 2 / 2 ). It plays a 
key role in the so called ‘dual cascade’ of two-dimensional turbulence, where nonlinearity 
generally results in a direct cascade of the spectral enstrophy density, enabling an inverse 
cascade of the spectral energy density (Fjprtoft 1953). This inverse cascade is often 


associated with the emergence of large-scale flow features, a general prediction of the 
MRS theory discussed in §1. 


Beyond the conservation of an infinite number of Casimirs, Eq. (2.1) further implies 


the existence of a topological constraint that contours of vorticity cannot cross. This 
nonlocal constraint is usually ignored in the construction of statistical descriptions. 


2.3. Numerical methods 

Two very different numerical methods are used to cross verify our conclusions. The first is 
a purely grid-based method using a quasi-isotropic spherical ‘geodesic’ grid consisting of 
D cells — hereafter referred to as the ‘geodesic grid method’ or ‘GGM’ (Heikes & Randall 


1995a|6 Qi fe Marston|2014 ). In the absence of explicit numerical dissipation, the GGM 
conserves both energy E and enstrophy Z. To remove enstrophy cascading to small scales, 
quad-harmo nic h yperviscosity is applied by adding the term z/ 4 (V 2 + 2 )V 6 cj to the right 
hand side of ( |2.l| ) (the (V 2 +2) operator ensures L is conserved). The coefficient z/ 4 for any 
given resolution D is chosen so that the most rapidly dissipating mode decays at a rate 
of 160 (initially, |cj| max = 5.8, see below). Adding explicit dissipation of course changes 
the original problem, but we argue it is necessary to account for the forward cascade of 
enstrophy from resolved to unresolved scales. Hyperviscosity, however, has unwanted side 
effects, and in particular it can amplify vorticity extrema ( Mariotti et ai [1994 


Jimenez 


1994 Yao et a/T||1995 ). Nonetheless, hyperviscosity is more scale selective than ordinary 
viscosity and much less dissipative overall for a given resolution D. Time stepping is 
carried out with a second-order leapfrog scheme, using a time step of At = 0.005 at 
all resolutions. (The time step respects the CFL stability constraint.) Decorrelation of 
even and odd time levels in the leapfrog scheme is avoided by applying a weak Robert- 
Asselin time filter (a = 0.001), improved as recommended in Williams] (2009), with the 
modification parameter 0.53. We have checked that using a smaller value of a = 10 -4 
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makes little difference. Hyperviscosity is implemented explicitly. GGM is freely available 
online jf] 

The second method is the ‘Combined Lagrangian Advection Method’ (CLAM, Dr itsche| 
|fc Fontane| ( [20 10)). This is a hybrid method, using material vorticity contours (solving 
dr/di = u{r,t )) to accurately and straightforwardly satisfy the conservative form of 
B but also making use of conventional grid-based procedures to invert \2.2\ and to 
interpolate u at contour nodes. CLAM is an extension of the closely related Contour- 
Advective Semi-Lagrangian algorithm (CASL, Dritschel & Ambaum (1997)). CLAM 
additionally solves (2.1) for uj on a grid and blends this grid-based solution with the 
contour-based one at the end of every time step to significantly improve energy conserva¬ 
tion compared to CASL. This is essential for carrying out exceptionally long simulations, 
as required here. In contrast to GGM which represents the vorticity as a continuous 
variable on a discrete grid in space, in CLAM the vorticity itself is discretized, and the 
contours can move continuously in space. A contour spacing of Ac o = 0.1 is used in all 
simulations here, corresponding to 109 contour levels. Details of the method along with 


tests and comparisons with other methods may be found in Dritschel & Fontane (2010). 


What is not described there is the specific implementation of CLAM in spherical geom¬ 
etry. Following Mohebalhojeh & Dritschel ( |2007 ), we employ a regular latitude-longitude 
grid of dimensions n g x 2 n g , enabling FFTs in longitude. Polar points are avoided by 
placing latitudes at — rc/2 + (j — l/2)n/n g , for j = 1, 2, ...n g . Latitude derivatives are 
computed spectrally after carrying out FFTs along great circles passing through the 


poles. The inversion of Poisson’s equation (2.2) for if is performed using 4th-order com¬ 


pact differencing in latitude for each longitudinal wavenumber. For the time evolution, 
CLAM evolves side-by-side a contour representation of cj, a gridded representation, and 
a residual gridded field to minimise numerical dissipation and to achieve excellent conser¬ 


vation properties (see Dritschel & Fontane (2010)). The gridded fields are evolved semi- 
spectrally, using the pseudo-spectral method whereby nonlinear products are carried out 
in physical space. Only the small residual field is damped by a weak bi-harmonic hyper¬ 
viscosity, with a damping rate of 2oj rms (t) on wavenumber n g (initially, uj rms = 1.7060324 
and decays thereafter). A 4th-order Runge-Kutta method is used for time stepping, with 
an adaptive time step At = min{0.l7t/|a;| ma x, 0.77r/(n^|i/| max )}. Contours are regularised 
periodically by ‘contour surgery’ at the subgrid scale n/ldn g (Dritschel 1988), and oc¬ 


casionally by contour regeneration, using the standardised parameter settings set out in 


Fontane & Dritschel (2009) 


Notably, the effective resolution of CLAM is approximately 16 times finer in each 
direction than the basic ‘inversion’ grid of dimensions 2 n g x n g , based on previous direct 
comparisons with standard pseudo-spectral methods ( Dritschel fc Scott||2009 
fc Fontane||2010 Dritschel fc Tobias||2012 ). 


Dritschel 


3. Results 


3.1. Flow initialisation and evolution 

A number of simulations differing only in resolution were carried out using both the GGM 
and the CLAM methods. In GGM, we used resolutions D = 163842, 655362 and 2621442, 
while in CLAM, we used resolutions n g = 128, 256, 512 and 1024. All simulations start 
with an identical vorticity field o;(r,0) built from a set of spherical harmonics with 
degree 4 ^ ^ 10; this field is shown in figure [l] (see appendix A for the amplitudes 


f GGM is implemented in the application “GCM” that is available for OS X 10.9 and higher 
on the Apple Mac App Store at URL http://appstore.com/mac/gcm 
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Figure 1. Initial vorticity distribution cj(r, 0) as seen from the equatorial plane (a) at 0° longi¬ 
tude, and (b) at 180° longitude. Colours (online) range from blue/purple for uj — c^min = —5.1187 
through green to red for uj = a; max = 5.7996. An orthographic projection is used. 



Figure 2. Early-time vorticity distribution cj(r, 4) shown as a function of longitude and latitude 
for (a) GGM at D m 655362, (b) GGM at D = 2621442, (c)CLAM at n g = 128, and (d) CLAM 
at rig = 256. The images are rendered as described in hgurefT] The percentage decay in enstrophy 
Z by this time is 0.13%, 0.02%, 0.48% and < 0.01%, respectively. Same color scale as in Fig. [l| 


used). Other initial conditions were tried and produce qualitatively similar results. The 
subsequent evolution, for early times, is nearly indistinguishable in all of the numerical 
simulations conducted — this is illustrated in figure [2] at t = 4 (note, the eddy turnover 
time 4ddy = 47r/|cj| max = 2.1668). Here, the two highest resolution simulations carried 
out using GGM are compared with the two lowest resolution simulations carried out 
using CLAM. By this stage of the evolution, the scale cascade of uj is still well resolved, 
though a small fraction of the enstrophy Z (< 1%) has already been lost. 























7 


On the late-time behaviour of a bounded, inviscid two-dimensional flow 
(a) (b) 



Figure 3. Early-intermediate time vorticity distribution cj(r,40), depicted as in figure 15] for 
(a) GGM at D = 655362, (b) GGM at D = 2621442, (c) CLAM at n g = 128, (d) GLAM 
at n g = 256, (e) CLAM at n g — 512 and (f) CLAM at n g — 1024. The percentage decay in 
enstrophy Z by this time is 56.8%, 52.2%, 50.0%, 44.1%, 38.3% and 32.4%, respectively. In 
CLAM, the results are plotted on a grid 16 times finer in each direction except in (e) where 
every other fine grid point is plotted, and in (f) where every 4th grid point is plotted. Same 
color scale as in Fig. [l] 


By t — 40, ten times further into the evolution, the simulations begin to diverge 
- see figure [3] which includes now the two highest resolution CLAM simulations at 
n g = 512 and 1024. By this time, significant enstrophy dissipation has occurred, albeit 
least in the highest resolution case. Vorticity contours are stretched nearly exponentially 
(at a rate comparable to |cj| max , see e.g. Dritschel et al. (2007)). By incompressibility, 
this forces filaments to thin exponentially and gradients to grow exponentially. This 
behaviour cannot be followed indefinitely in numerical simulations at finite resolution, so 
inevitably dissipation occurs. In CLAM, unlike in GGM and other grid-based numerical 
methods, the steep gradients can be preserved, without dissipation, due to the contour 
representation of uo used in CLAM. This preservation of steep gradients is the key feature 
which enables CLAM to model flows at very long times albeit at the cost of representing 
the continuous vorticity field by a discrete set of levels. There is no progressive erosion of 
vorticity gradients, as in grid-based methods. The vorticity filaments, on the other hand, 
play a relatively benign role in the dynamics. They contribute negligibly to the flow field 
and have little chance of rolling up into coherent vortices (Dritsch el et a/.||1991 ). 
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(a) (b) 



Figure 4. ID vorticity power spectra Z(l,t) at (a) at t — 4 and (b) t — 40. Colours (online) 
cycle from black, to blue, to red, to green and back to black, starting from the lowest resolution 
GGM simulation (D = 163842) and ending with the highest resolution CLAM one (n g = 1024). 
Three spectra are shown for GGM and four are shown for CLAM. Log (base 10) scales are used. 


Notably, the CLAM simulations exhibit more structure than the GGM ones. The 
GGM simulations appear to be at lower resolution, and the effects of numerical diffu¬ 
sion are visible. The ID vorticity power spectra Z(l,t), shown in figure [4] at time t = 4 
and 40, confirm this impression. The CLAM spectra are undissipated out to the maxi¬ 
mum wavenumber computed (which is limited to spherical wavenumber £ = 2800 in the 
SHTOOLS software package used). The GGM spectra turn down from hyperdiffusion 
but also exhibit a rise near the maximum wavenumber (the spectral amplitudes are how¬ 
ever too weak to play a significant role). Note the significant shallowing of the spectrum 
in time, as nonlinearity transfers a large fraction of the enstrophy to small scales. The 
spectral slope over the range Z ^ 100 shallows from approximately —2.9 to —1.9 between 
t — 4 and 40 in the highest resolution CLAM simulations. 

By t — 400, a further ten times later in the evolution, much of the complex structure 
evident at t = 40 has cascaded beyond the maximum scale resolved in all simulations 
— see figure [5] for the two highest resolution GGM and CLAM results. Nonetheless, 
extensive small-scale structure survives, particularly in the CLAM results. The main 
quadrupolar vortices exhibit stepped distributions, ‘vorticity staircases’ (evident also in 
the GGM results), together with a multitude of small-scale vortices caught between the 
sharp vorticity gradients. This structure persists, albeit somewhat reduced, all the way 
to the end of the simulations at t — 4000 (see figure [6]). A zoom at t = 4000, shown 
in figure [7] for the highest resolution CLAM case, reveals the extent of this fine-scale 
structure — literally hundreds of small-scale vortices exist between the vorticity steps 
on this large-scale vortex alone. Notably, virtually all of these small-scale vortices have 
vorticity anomalies which cause the vortices to rotate with the shear associated with 
the differential rotation of the main large-scale vortex they are embedded in. This gives 
rise to stabilising ‘cooperative shear’, enabling them to persist for long times, perhaps 
indefinitely ( Dritschel] 1990). 

The reduction in the size and number of small-scale vortices between t = 400 and t = 
4000 is likely to be a numerical artefact, the effect of very weak dissipation accumulated 
over a very long time. In CLAM, every few eddy turnover times (about 7 or so units 


of time), contours are rebuilt on an ultra-fine grid to prevent overlapping (Fontane & 
|PritscheT |2009 Dritschel & Fontane 2010| ). Each time this happens, inevitably small 
contours tend to shrink slightly on average (and the effect is strongest for the smallest 
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Figure 5. Intermediate-late time vorticity distribution ct;(r,400), depicted as in figure [ 2 ] & [ 3 ] 
for (a) GGM at D = 655362, (b) GGM at D = 2621442, (c) CLAM at n g = 512 and (d) GLAM 
at n g — 1024. Note, only the two highest CLAM results are shown, but the lower resolution 
results are qualitatively similar. By this time, the enstrophy Z has decayed to approximately 
21% of its initial value in all cases. No further significant decay occurs to t — 4000 (< 0.05%). 
Same color scale as in Fig. [I] 


(a) (b) 



Figure 6. Late (and final) time vorticity distribution ^(r, 4000), depicted as in figure [ 2 ] h [ 3 ] 
for (a) GGM at D = 655362, (b) GGM at D = 2621442, (c) CLAM at n g = 512 and (d) GLAM 
at n g = 1024. Same color scale as in Fig. [l] 
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Figure 7. Zoomed, full-resolution view at t — 4000 of the vortex seen in the lower-left part of 
figure[6jd). The view is centred at longitude —157r/16 and latitude — 7r/8, and shows a 7r/2 range 
in longitudes and latitudes (l/8th of the domain or 8192 2 grid points). Colours (online) range 
from the minimum to the maximum vorticity in the domain of view. Concentric steps seen in 
the vorticity are an artifact of the discretization of the vorticity levels in the CLAM algorithm. 


resolved contours). This, we believe, is the primary cause for the reduction of small-scale 
structure. Notably, each resolution doubling performed here reveals significantly more 
structure persisting to the end of the simulation, and in the limit of infinite resolution 
- a perfect model — we expect a multitude of vortices occupying an ever extending 
range of scales. Unlike in the case of inviscid 2D turbulence in an unbounded space 
(Dritschel et al. 2009), where vortices are free to roam and occasionally collide to form 
both larger and smaller vortices, here the large-scale flow imposed by the finite domain 
(the sphere) channels vortices to move along the quasi-steady streamlines of each large- 
scale vortex, preventing vortex collisions. As a result, vortices may persist indefinitely 
once the possibility of encountering another vortex in a nearby streamline at close range 
has been eliminated. 

Vorticity cross-sections through the 4 main vortices in figure |6|d) exhibit stepped vor- 
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(a) (b) 




Figure 8. (a)-(d) Constant-latitude cross sections of the vorticity field taken through each of 
the 4 large-scale vortices (from left to right in figure [6jd) ) in the highest resolution CLAM 
simulation at t — 4000. The relative longitudinal grid point is indicated on the horizontal axis. 
The small steps show the vorticity contour level spacing used in the CLAM simulation. The 
important features are the large steps, seen in many of the vortices. Note also that the cross 
sections occasionally pass through a small-scale vortex. 


ticity distributions, with alternating steep and shallow gradients (see figure [8]). These are 
also found in the highest resolution GGM simulation, together with small-scale vortices 
(though many fewer due to the effects of hyperdiffusion). The stepped vorticity distribu¬ 
tions seen here appear to be closely analogous to those arising spontaneously in rotating 


2D turbulent flows having a mean ‘planetary’ vorticity gradient (see Dritschel & McIn¬ 


tyre (2008) for a review). There, turbulence cannot entirely break down this gradient 


but rather mixes vorticity inhomogenously, resulting in a stepped vorticity distribution 
in the inviscid limit (Scott & Dr itschel| [2012). The associated flow consists of a series of 
eastward (prograde) jets centred on each step, reminiscent of the jets observed in many 
planetary atmospheres, most notably Jupiter’s. In the simulations conducted here, each 
large-scale vortex has a vorticity gradient which is reshaped into a stepped distribution 
early on in the turbulent flow evolution. 

The impact of this turbulent mixing can also be seen in the area fraction occupied 
by each vorticity level — the so-called ‘vorticity measure’ p(uj,t). This is shown at a 
few selected times in figure [9] for the highest resolution GGM and CLAM simulations, 
in (a) and (b) respectively. In a perfectly inviscid fluid, p(cj,£) is time-independent and 















12 

(a) 


D. G. Dritschel, W. Qi and J. B. Marston 
(b) 




Figure 9. The log-scaled vorticity measure log 10 p(cj, t) at selected times in (a) the highest 
resolution GGM simulation, and (b) the highest resolution CLAM simulation. The times shown 
(colours online) are t — 0 (black), t = 40 (blue), t = 400 (red) and t — 4000 (green). 


generates the conserved Casimirs in (2.5) through the relations C n = 4n f cj n p(u;)du;, 
n = 2, 3, .... The lack of conservation seen here is caused by the inevitable cascade of 
vorticity to small scales — none of the Casimirs can be preserved in simulations at finite 
resolution (indeed, C 2 decays by 79% here). The initial vorticity measurep(cj, 0) is broadly 
distributed between u; m j n = —5.1187 and u: max = 5.7996, then narrows and diminishes 
nearly everywhere except for values of uj near 0. By t = 400, the vorticity measure has 
converged to a nearly fixed form, which is noticeably jagged, perhaps related to the 
non-monotonic vorticity profiles seen in figure |8ja). In the CLAM results on the right, 
notice that cj m j n and co max are conserved, and that the area occupied by these extreme 
vorticity values does not change significantly over time. This is sensible physically, as the 
extreme values are most resilient to deformation. On the other hand, in the GGM results, 
there is clear evidence of hyperviscous overshoot: the extreme values are not conserved 
(at times not shown here they are nearly twice their initial values). Slightly greater and 
unphysical extreme values persist to the end of the simulation, though they have very 
small measure. In both simulations, the observed strong growth of p(cj,£) near uj = 0 
is due to the efficient mixing of weak vorticity, which essentially behaves as a passive 
scalar. As \co\ increases, it becomes increasingly difficult to shear out vorticity structures 
and mixing becomes less efficient. 

3.2. Unsteadiness 

We next address the possible long-term equilibration of the flow, a core prediction of 
the MRS statistical-mechanical theory. As t 00 , that theory predicts a steady flow 
characterised by a functional relation between vorticity and streamfunction u) = F('iji). 
The precise functional F depends on the energy, angular momentum (here zero), cir¬ 
culation (here zero due to the absence of boundary) and the unobservable higher-order 
fine-grained Casimirs. The resolved ‘coarse-grained’ Casimirs associated with the non¬ 
cascading vorticity field are often used to approximately determine the values of the 


fine-grained Casimirs (see e.g. Qi & Marston (2014) & references therein). The details do 
not matter for what concerns us here. If the flow is steady, then F can be obtained from a 
scatter-plot of ^ vs cj, in which F would appear as a single curve. Our results are shown 
for both GGM and CLAM at the highest resolutions and at the final time in figure [lO] 
First of all, there are clearly multiple branches, each corresponding to one of the main 
vortices. Second, the data do not collapse on a line, but rather in a band, indicating that 
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(a) (b) 






if 


Figure 10. Scatter plots of vorticity w vs. streamfunction at t = 4000 in (a) the highest 
resolution GGM simulation, and (b) the highest resolution CLAM simulation. 


(a) (b) 




Figure 11. Early time evolution (0 ^ t ^ 80) of (a) the ‘configuration’ W{t) (defined in the 
text) and (b) the rate of enstrophy dissipation — dZ/dt, for both the highest-resolution CLAM 
and GGM simulations. Colours (online) show CLAM in black and GGM in blue. 


the flow is not steady (see below for further evidence). Third, in the CLAM results, one 
sees a myriad of small-scale vortices, due to the much weaker level of dissipation present. 
In fact, due to the need to keep the figure size manageable, the CLAM results are down- 
sampled by a factor of 16 in each direction. The true number of small-scale vortices is far 
greater. Figure (lO^b) gives a greatly simplified view. The upshot is that the flow remains 
unsteady and extremely complex, even despite the inevitable dissipation in CLAM which 
acts to remove small-scale structure over time (see discussion above). 

An alternative way of measuring the unsteadiness of the flow was introduced by |Qi fc| 


Marston (2014). They developed a rotation and amplitude independent measure of the 
flow configuration that involves only the streamfunction amplitudes of the spheri¬ 
cal harmonics Y™ of degree i = 2 (these coefficients characterise the main quadrupolar 
pattern seen at late times). The order m satisfies — 2 < m < 2. Using the standard con¬ 
vention for spherical harmonics (for which = (—l) m y^ m * where * denotes complex 
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Figure 12. Late time evolution (t ^ 200) of the ‘configuration’ W(t) for both the second 
highest resolution GGM simulation (top) and highest resolution CLAM simulation (bottom). 


conjugation), the configuration is obtained from the amplitudes: 


W(t) = 


-2 V 4 + 6020 ( 02,-1 021 + 2 02,-2 0 : 22 ) ~ 3^6 (' 02-2 021 + 02,-1 022 ) 
( 02,-2 + 02,-1 + 020 + 021 + 022 ) 2 


(3.1) 

W ranges from —2 to +2; when \ W\ = 2, there are just two vortices, while when W = 0 
there are four of equal intensity. The sign of W indicates the sign of the dominant pair 
of same-signed vortices. 

The evolution W(t) at early times, 0 ^ t ^ 80, in both GGM and CLAM at the 
highest resolutions is compared in figure [llja). The correspondence is excellent until 
around t = 20, a little after the time of peak enstrophy dissipation rate — dZ/dt (shown 
in (b)). After this time, the results diverge due to the loss of predictability. However, 
both methods — at all resolutions — exhibit the same qualitative behaviour in W(t) 
thereafter: quasi-periodic and non-diminishing oscillations, i.e. persistent unsteadiness. 
This is shown in figure [T2| for two very long integrations carried out to around t = 10000; 
the top panel shows GGM at D = 655362 while the bottom shown CLAM at n g = 1024 
(both for t ^ 200). Movies of the flow evolution indicate that the oscillations are caused 
by the four main vortices moving around each other, with global-scale excursions, and 
with no sign of relaxing to equilibrium. Note that similar oscillations are also seen in the 
dipole case of figure 2 of Morita (2011) where they are described as ‘small perturbations’ 
without detailed discussion. 

A simple model of this vortex motion is to consider just 4 point vortices, singular dis¬ 
tributions, in place of the continuous vorticity distribution above. For point vortices, the 
dynamical system has just 3 degrees of freedom, after taking into account conservation of 
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Figure 13. Final 10 periods of evolution (approximately 192 time units) in the CLAM, GGM 
and point-vortex simulations. Time is shown relative to t max — 192 in each simulation. Colours 
(online) show CLAM in black, GGM in blue, and the point vortex model in red. 


energy, angular momentum and the fact that the dynamics depends only on the relative 
(angular) separation of the vortices (Polvani & Dritschel 1993; Newton|[2001 ). The vor¬ 
tex motion is therefore either quasi-periodic or chaotic in general; steady configurations 
are rare (indeed of zero probability for arbitrary choices of the vortex circulations and 
positions). 

To test the effectiveness of this simple model, the 4 main vortices from the highest 
resolution CLAM simulation at t = 4000 (see figures [6] and [8| were identified by locating 
contiguous vorticity regions having \oj\ > 0.2, a small threshold value enabling one to 
capture most of the circulation within each vortex. The circulations (k = 1, 2, 3 & 
4) were found by area integration over each contiguous region of vorticity Rk, and the 
centres from 


r k = f [ [ ur dfi. (3.2) 

1 k J jR k 

This procedure however does not lead to zero angular impulse L , which for point vortices 
is equal to k Since L = 0 in the original flow, we adjust the circulations T^ by 

subtracting A • rk for each k, where A is determined from a 3 x 3 linear system. This 
adjustment is not biased toward any one vortex and leads to a 10 to 15% change in the 
circulations. 

With Tk and rk determined in this way, the governing dynamical system 


drk _ 1 
d t 4tt 


47r 




Vj x r k 

1 - r j ■ r k 


(k = 1, 2, 3, 4) 


(3-3) 


(P olvani fe Dritschel| l993) was integrated forwards for 4000 time units, with the results 
diagnosed as above for the configuration W(t) to compare with the CLAM and GGM 
results. The key finding is that W(t) exhibits the same behaviour, in particular with 
the same dominant frequency of oscillation; the last 10 periods of oscillation in all three 
simulations is shown in figure [l3j (The averages and amplitudes of W(t) differ slightly 
due to the largely unpredictable nature of turbulence and the crude fitting method of 
the point-vortex system.) Much longer point-vortex simulations carried out to t = 10 6 
exhibit the same behaviour, with undiminished oscillations. 

The correspondence between the three simulations is even more remarkable if one 
examines the frequency power spectrum of W, namely W(f) = \W \ 2 , where W is the 
Fourier transform of W and / is the frequency (here equal to the inverse period). The 
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Figure 14. Configuration frequency spectrum VV(/) for the second highest resolution GGM 
simulation (top), the highest resolution CLAM simulation (middle), and the point vortices (bot¬ 
tom). In GGM, W(/) was computed from the time series of W for t G [400,10000]; in CLAM 
we used [400,9022], while in the point-vortex model we used [0,4000]. Note: log 10 W is plotted 
versus log 10 /. 


frequency spectra (log 10 scaled) for the three simulations are shown in figure [l4j vertically 
aligned, with GGM on top, CLAM in the middle, and the point vortices on the bottom. 
Most strikingly, the dominant frequency f — 1/T with T ~ 19.14 is matched in all 
three simulations. The power at this frequency is orders of magnitude above surrounding 
frequencies, indicating that the dynamics are very nearly periodic. Other less prominent 
peaks also line up, corresponding to longer period oscillations. Overall, the agreement is 
remarkable, given the enormous simplifications made to reduce the dynamics to that of 
4 point vortices. Most importantly, the point vortex model strongly supports our claim 
that the dynamics never comes to equilibrium. 


4. Conclusions 

In this paper, we have studied the long-time asymptotic behaviour of a bounded two- 


dimensional inviscid flow. Statistical-mechanical theories stemming from Miller (1990); 


Miller et al. (1992) and Robert & Sommeria (1991); Robert (1991) argue that this be¬ 


haviour reaches a steady equilibrium state, characterised by a functional relation between 
vorticity and streamfunction. We have tested these theories by carrying out very long 
time simulations of 2D turbulence on a sphere, for zero angular momentum, employing 
highly distinct numerical methods. The key finding, which is robust across both methods 
at all resolutions, is that the flow remains persistently unsteady with no indication of 
equilibration in a steady flow. 

Of course, we cannot carry out simulations indefinitely, and one can always argue that 
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equilibration will happen eventually. However, we have shown there is a plausible alter¬ 
native. At long times, the vorticity field is primarily concentrated in 4 large vortices. 
They are generally unequal in circulation, and not regularly spaced. If one replaces each 
vortex by a point vortex of the same circulation, we have shown that the point vortices 
typically exhibit unsteady quasi-periodic motion — persistent unsteadiness. Fundamen¬ 
tally, we argue, the observed unsteadiness in the Euler equations for continuous vorticity 
is intimately associated with the unsteadiness of 4 point vortices moving on the surface 
of a sphere. 

We have also found that small-scale vortices pepper the flow field at high resolution. 
These vortices commonly have vorticity anomalies which cause the vortices to rotate 
in the same sense as the shear induced by differential rotation of the main vortex they 
circulate about. The small-scale vortices therefore experience a stabilising cooperative 
shear (Dritschel 1990). As a result, even these small-scale vortices may persist indefinitely. 
Moreover, as the resolution increases, more structure at small scales appear. Our main 
conclusion is that bounded inviscid flows on a sphere, with zero angular momentum, 
typically do not settle down to an equilibrium; instead, a few robust compact vortices 
emerge from the early flow evolution and interact chaotically or quasi-periodically. 

In the much studied case of doubly-periodic geometry, by contrast, often two opposite- 
signed vortices emerge. In this case, the corresponding point-vortex system is integrable, 
and the only possible motion is simple steady translation. This difference from spherical 
geometry greatly facilitates the emergence of a steady flow (in a translating frame of 
reference), though one still cannot rule out the persistence of small-scale vortices circu¬ 
lating around the main large-scale vortices. Inviscid flows in a doubly-periodic domain 
are likely to be less unsteady than flows on a sphere, though unsteadiness may persist 
forever. 
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KITP is supported in part by the NSF Grant No. NSF PHY11-25915. This work was 
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Appendix A. Initial Condition 

The initial streamfunction at time t = 0 is a superposition of spherical harmonics with 
spherical wavenumber in the range 4 ^ i < 10, 

10 i 

m ^)=EE Y t m (e, 0), (a i) 

1=4 m=—t 

with ipe-m = ( —l) m, 0| TO and the following choice of amplitudes ipe m : 

v >40 = -0.0114727 
V>4i = 0.0255012 + 0.0337679 i 
v> 42 = -0.0366116 - 0.0312419 i 
ip 43 = 0.000614585 + 0.000602384 i 
V>44 = 0.0160784 - 0.0194655 i 
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0 5O = 0.00357198 
0 51 = -0.000446705 - 0.0100286 i 
0 52 = 0.0100984 + 0.00778139 i 
0 53 = -0.0142904 - 0.00864324 i 
054 = 0.0156104 - 0.0111872 i 
-055 = -0.00275238 - 0.0106832 i 

ipso = -0.00352493 
0 61 = -0.00739953 + 0.00218422 i 
062 = 0.0140976 + 0.00864871 i 
0 63 = -0.0141289 - 0.00654953 i 
064 = 0.0200989 - 0.0152088 i 
065 = 0.00682398 + 0.0132045 i 
066 = 0.00919374 + 0.00889445 i 

0 7 o = 0.00269739 
0 71 = -0.00607359 - 0.00179029 i 
072 = 0.00299202 - 0.0115332 i 
0 7 3 = -0.00418447 - 0.000665072 i 
074 = 0.00853673 - 0.0011805 i 
0 75 = -0.00229449 - 0.00731682 i 
076 = 0.00315312 - 0.00332254 i 
077 = 0.0127047 + 0.000153664 i 

ipso — —0.0070908 
0 81 = -0.00820113 - 0.00180399 i 
0 82 = 0.00726466 - 0.00109238 i 
083 = 0.00327615 + 0.00732665 i 
0 84 = -0.00418371 + 0.00764604 i 
0 85 = -0.00506492 + 0.00152613 i 
0 86 = 0.00289222 - 0.0126644 i 
0 87 = -0.00113653 - 0.00983086 i 
088 = 0.00390903 + 0.00337946 i 

090 = 0.00366743 
09! = -0.00670701 - 0.0052387 i 
0 92 = -0.00102863 - 0.00831521 i 
0 93 = -0.00429645 + 0.00881862 i 
094 = 0.00487287 - 0.00708109 i 
0 95 = -0.00183309 - 0.00476393 i 
096 = 0.0124893 - 0.00170479 i 
097 = 0.00438044 + 0.00091529 i 
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^98 = -0.00309726 - 0.000427295 i 
^99 = -0.0053439 + 0.00203709 i 

Vhoo = -0.00507863 
Vhoi = -0.00247735 + 0.00024564 i 
V402 = -0.00357194 - 0.00141963 i 
Vhos = -0.00077508 + 0.00719427 i 
Vho4 = -0.00462811 - 0.00550109 i 
Vhos = -0.00101044 - 0.00773111 i 
Vhoe = 0.00395564 + 0.0027653 i 
Vho? = 0.00121852 + 0.00310834 i 
Vhos = 0.00314655 + 0.000289766 i 
^io9 = -0.00279021 - 0.00314691 i 

V’ioio = 0.00258861 - 0.000254598 i (A 2) 
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